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Abstract 



^^ ^he wall shear stress is a quantity of profound importance for clinical diagnosis of artery diseases. The lattice Boltzmann is an easily 
H parallelizable numerical method of solving the flow problems, but it suffers from errors of the velocity field near the boundaries 
^^ 'which leads to errors in the wall shear stress and normal vectors computed from the velocity. In this work we present a simple 
<^ formula to calculate the wall shear stress in the lattice Boltzmann model and propose to compute wall normals, which are necessary 
T-l- ito compute the wall shear stress, by taking the weighted mean over boundary facets lying in a vicinity of a wall element. We carry 
T— I out several tests and observe an increase of accuracy of computed normal vectors over other methods in two and three dimensions. 

XJsing the scheme we compute the wall shear stress in an inclined and bent channel fluid flow and show a minor influence of the 

(— I normal on the numerical error, implying that that the main error arises due to a corrupted velocity field near the staircase boundary. 
Q_iFinally, we calculate the wall shear stress in the human abdominal aorta in steady conditions using our method and compare the 



Oh 



results with a standard finite volume solver and experimental data available in the literature. Applications of our ideas in a simplified 
protocol for data preprocessing in medical applications are discussed. 

fluid flow, numerical methods, the lattice Boltzmann method, LBM, wall shear stress, WSS 
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Computer simulations of the blood flow in the human circu- 
lar system allow physicians to visualize the internal flow struc- 
ture, a capability that might have a significant impact on diag- 
nosis and medical treatment of arterial diseases (for review see 
jlll|2l|3|]). Digital simulations have been applied in many sce- 
narios and for various purposes e.g. to study the role of hemo- 
dynamic forces in the development of atheromatous plaques in 
human cartoids atherosclerosis ||4|, in a clinical study on the 
rupture risk assesment in the cerebral aneurysm flSl or in a dis- 
covery of a shear-driven mechanism of platelet aggregation in 
arterial thrombosis |6i|. Direct results of hydrodynamical simu- 
lations (velocity and pressure fields) are not practical in medi- 
cal analysis. To better understand the impact of the blood flow 
dynamics on the complex biochemical phenomena associated 
with the development of vascular diseases, a wall shear stress 
(WSS, see Fig. [TJ with its gradient (WSSG), a hoop stress or 
an oscillatory shear index (OSI) are often computed. WSS is a 
tangential force exerted on the unit of the wall surface by the 
flow im and OSI describes its deviation from the average direc- 
tion (si] . WSS has been a subject of intensive research for many 
years, includingstudies on its role in development of human 
atherosclerosis ||9lll0|], the growth of an intracranial aneurysm 
I n|] and development of the dissecting aneurysm in the human 
aorta 11211 . WSS determines the structure, function and gene 
expression of an endothelial cells ||4J, ll3[]. It is now well es- 
tablished that a low shear stress is correlated with rapid devel- 
opment of various vascular malfunctions lll4l llSi Il6i Il7i llSll . 
These facts motivate physicians to find applications of the WSS 




Figure 1: (color online) Wall shear stress magnitude in the bifurcation of the 
human abdominal aorta (the lighter color the higher the stress) calculated using 
a lattice Boltzmann model (left) and the finite volume method (right). 



in diagnosis 0191 12011 and to seek efficient and practical proto- 
cols for its measurement. 

In order to compute WSS for patient-specific data one has to 
acquire the organ geometry, convert and import it to the Navier- 
Stokes solver, compute the velocity field and take gradients of 
its components at the wall. The process may be iterative and 
further time-averaging might be necessary if the flow is tran- 
sient. The major difficulty of this procedure is the prohibitively 
long pre-processing and computation time which limits its ap- 
plicability in a real patient diagnosis. The pre-processing in- 
volves mainly the conversion of digital data from the computer 
tomography scans into a three-dimensional mesh required by 
standard finite element/finite volume solvers 11211 l22l 12311 . High 
quality of the output mesh is required by flow solvers to work 
properly and computer-human interaction during the meshing 
process is often necessary. The flow solver runtime depends on 
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many factors e.g. solver type, mesh resolution, flow properties 
and hardware used. It may be a bottleneck of the whole proce- 
dure as well, particularly for geometries of high resolution, time 
dependent flows and problems involving fluid-structure interac- 
tion. It is widely believed that the run time of fluid solvers could 
be significantly reduced by using efficient parallel algorithms 
on massively parallel architectures, e.g. graphics processing 
units (GPUs) m\ . 

Time issues mentioned above might be addressed with the 
lattice Boltzmann method (LBM), a recently developed method 
for solving fluid flow problems which is based on the kinetic 
theory of gases 1251]. The method operates directly on voxel 
data and thus does not require meshing at the level of accuracy 
required by standard solvers. It also proved excellent parallel 
scalability on supercomputers 11261 12711 as well as on emerging 
parallel architectures ll28[ |29[ BOi Bill . The method was vali- 
dated against analytical solutions in a two-dimensional oscilla- 
tory channel flow and its second order accuracy in space and 
time was confirmed l32l 13311 . In medical applications, LBM 
was tested in solving time-dependend three-dimensional flows 
in models of the human abdominal aorta ll34ll and compared 
against finite element solvers in the superior mesentric artery 
blood flow 13511 . Non-Newtonian LBM models also exist e.g. 
the Carreau-Yasuda fluid model applied in cerebral aneurysm 
problem 13611 . Palabos, a parallel open source LBM implemen- 
tation, was used to model blood flow in cerebral aneurysm and 
compared to a finite volume flow solver OTJ 13811 and the differ- 
ence between the velocity fields obtained in both methods was 
estimated to be around 10%. AppUcation of LBM to modeling 
aneurysm and tumor development were recently discussed 113 911 . 
He et al. compared LBM combined with the level-set method 
to an OpenFOAM solver in patient-specific geometry ll40ll . 

Despite the second order accuracy of the LBM stress ten- 
sor in bulk fluid being recently reported ll4l[ 14211 . the WSS in 
LBM suffers from inaccuracy introduced by a staircase bound- 
ary. The problem may be eased with subgrid boundary con- 
ditions 14311 or the recently developed unstructured mesh LBM 
formulation ll44ll . but both of them introduce additional com- 
plexity into the solver code, making it more difficult to main- 
tain and parallelize. Recently Stahl et al. showed that a few 
lattice nodes away from the wall the error is significantly de- 
creased in the standard LBM formulation ll45ll . He suggests to 
calculate normal vectors required by this procedure from the 
velocity field as the fluid flowing along the noslip walls follows 
its geometry. Thus, the normal at the boundary is approximated 
by the normal to the local velocity vector in a cell close to the 
wall. However, this method behaves peculiarly bad in a vicin- 
ity of the wall, with maximum error close to 90%. Moreover, 
it is not suitable for time dependent flow problems unless one 
solves an additional steady flow case. It also requires a solution 
to the eigenvalue problem and does not provide the vector sense 
relative to the wall. 

In this paper we discuss another way of computing the nor- 
mal vectors which does not suffer from above issues. If we 
look at the LBM staircase wall, we intuitively see its orien- 
tation. This is because we subconsciously average the input 
of all faces we see. A similar example is our brain's han- 



dling of a two-dimensional line on a computer display: if pixels 
are small enough, the line orientation is easily recognized by 
eye. We propose to utilize similar averaging of a few boundary 
facets around the point of interest to compute the normal vec- 
tors required for evaluating WSS. The procedure is simple and 
its implementation consists of only a few algebraic operations. 
Therefore, no eigenvalue problem has to be solved. Addition- 
ally, the information about the normal sense is provided, and 
we show that our procedure tends to provide more accurate re- 
sults. Henceforth the normals calculated in this way will be 
called 'geometric normals' and the normals approximated from 
the velocity field will be called 'dynamic normals'. 

The structure of the paper is as follows. In Sec. |2]and|3]we 
briefly introduce the lattice Boltzmann method and the formula 
for the wall shear stress. In Sec. I3.1l we explain how to com- 
pute geometric normals and verify them in simple geometries. 
The wall shear stresses in an inclined and a bent channel flow 
are calculated in Sec. |4] A comparison of WSS obtained with 
the lattice Boltzmann and a finite volume methods in the hu- 
man abdominal aorta is given in Sec.|5] Section |6]presents the 
discussion of the results and Sec. |7]is devoted to conclusions. 

2. The lattice Boltzmann method 

LBM is a computational method of solving the Navier- 
Stokes equations. It was developed mainly as a remedy to noise 
issues of the lattice gas automata (LGA) 143, |43]. In LGA the 
local state of the system is described by a boolean variable Ci, 
where i runs over lattice vectors (vectors that link neighbour 
sites in the lattice). In contrast, LBM replaces c\ by a parti- 
cle distribution function /i e [0, 1], which is interpreted as the 
probability that a particle found in the lattice node moves along 
the /-th direction. A typical LBM algorithm runs two steps, 
propagation and collision, and is expressed by a set of discrete 
equations: 



/„(x + Cadt, t + St)= fa(x, t) + Q„(f), 



(1) 



where c^^ is the a-th lattice vector, 6t is the time interval (usually 
dt = 1) and Q„(f) is a collision operator The collision operator 
may realize either a linear relaxation to the equilibrium distri- 
bution function or the multirelaxation of separate hydrodynam- 
ical and kinetic moments 14811 . The Chapman-Enskog analysis 
cane be used to show that the model reproduces incompressible 
Navier-Stokes equations [49(1. 

As an implementation of LBM we used the Sailfish library 
15011 . Sailfish works on GPU processors. It is written in Python 
and utilizes template-based methodology for automatic device 
code generation on CUDA and OpenCL GPU platforms. The 
original code is freely available under the LGPL licence. 

3. The wall shear stress 

When fluid flows over a rigid surface, the velocity at the wall 
vanishes and the no slip boundary condition holds. However, in 
a vicinity of the wall the tangential component of the velocity 
does not vanish; the corresponding gradient of the tangential 
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Figure 2: (color online) (left) A geometric normal vector (long arrow denoted 
as n) is calculated by averaging input from normals of neighbouring fluid/no- 
slip boundary nodes (short aiTows). Right: normal vectors computed for each 
individual boundary node. Geometric nomials are enlai'ged for easier viewing. 
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Figure 3: (color online) Geometry of the ball for L=25 (lattice units). Normal 
vectors calculated at half vertical distance are computed with Eq. for y = 1 . 
The averaging radius is r = 4 l.u. 



velocity component along the wall normal generates the wall 
shear stress, a force that is exerted by the fluid on the wall's 
surface. WSS may be expressed as the difference between the 
Cauchy stress on a plane and its projection on the plane normal 



(see Appendix A i. For an incompressible, Newtonian fluid the 
WSS in the lattice Boltzmann method can be calculated in terms 
of the non-equilibrium distribution function f^"'': 
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(2) 



where c, is the lattice sound speed, fi is the dynamic viscosity, 
g is the fluid density, and co is the LBM relaxation parameter. 
Ca.x and rix denote lattice and normal vector components, i.e. 
«, is the i-th component of the wall normal vector n and c„, is 
the i-th component of the lattice vector c„. Einstein summation 
convention is implied. A detailed derivation of Eq. [2] is given 
in [Appendix A| with its algorithmic form in [Appendix B| 

3.1. Wall orientation 

To calculate WSS with Eq. (O the wall orientation and the 
corresponding wall normal n must be known. A typical three- 
dimensional LBM simulation works on a discrete, uniform grid 
built of adjacent cubic cells. If a fluid node has a no-slip bound- 
ary neighbour then there is a boundary facet between them at 
which the no-slip condition of zero velocity is fulfilled. The 
facet is perpendicular to the line that joins node centers whereas 
its physical location depends on the LBM relaxation time. For 
a three-dimensional model there are six possible boundary nor- 
mals: +ex = (+1,0,0), +ej, = (0,+l,0), and+Cz = (0,0,+l). 

In order to calculate n we suggest taking an average of nor- 
mal vectors of the individual fluid/no-slip boundary facets in 
a spherical vicinity of the cell. For this we find all boundary 
facets whose centers lay in the sphere of radius r centered at 
a given cell (Fig. |2^). If b,- is the boundary normal of the j-th 
facet then the above procedure may be expressed as the arith- 
metic mean: 



„i,-i 



(3) 



3.2. Weighted scheme 

In the application of Eq. (O an unacceptable averaging of 
normals that belong to separate or opposite fragments of the 
object's surface may take place if r is taken to be too large. In 
the hemodynamical scenario, this may occur in small arteries 
or bifurcation regions where object radii are smaller that r and 
details may be lost due to averaging normals to opposite facets. 
Moreover, if r is too small, n is influenced by the staircase ap- 
proximation of the surface and may change rapidly from node 
to node. To deal with these issues we suggest to weight the 
mean in Eq. (O: 



n = ^ w, ^ w,b,- 



where nt is the number of facets within the neighborhood. The 
same procedure is repeated for each boundary cell (Fig. right) 



(4) 



where the sum goes over neighour facets. We choose w,- = 
1/(1 + diY, where c/, is the distance between the center of the 
node for which the average is computed and the center of the 
node with normal b,. 

The exponent y defines how strongly concentrated the 
weighting is, with y = corresponding to uniform weights (re- 
ducing to Eq. (|3]l); the strength of weighting increases with the 
value of y. 

3.3. Numerical tests for normal vectors 

We test the above procedure in three-dimensional ball and 
tube geometries for which exact normal vectors are known. 
First, a boolean array of size L^ is generated such that it holds 
if the voxel (three dimensional pixel) center lays inside the 
object and 1 otherwise. Boundary cells are represented by vox- 
els equal to 1 having at least one adjacent zero neighbour The 
object (ball or tube) radius is, R - L/2 and we test three res- 
olutions: L - 25^,50'' and lOO"*. We calculate normals on the 
circle created by boundary nodes at half vertical position. As an 
example, in Fig. |3]the geometry of a ball for L=25 is visualized 
with normal vectors calculated using our procedure. Fig. 0de- 
picts the maximum error between geomeric and exact normals. 
Results are compared to dynamic normals at various distances 
from the wall |45]. Geometric normals are found to be more 
accurate. To test the influence of the averaging radius on the 
accuracy of the normal for various exponents y we vary r and 
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Figure 4: Tlie maximum eiTor of a nornial vector (measured as the angle 
between numerical and exact vectors) for two numerical methods in a three- 
dimensional ball. Arrows represent geometric normals computed with Eq. J4) 
using y = 1 and r = 4. 
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Figure 5: (color online) The maximum eiTor of geometric nonnal J4) for three 
values of 7 as a function of the averaging radius r in a three-dimensionall ball 
and tube (inset) with radii R=25. 



calculate the error on a ball and a tube of radii R. The result- 
ing error is a non-monotonous function of r with a minimum at 
at r/R = 1 (see Fig. |5]l. At r/R > 1 the mean error for y = 
grows rapidly whereas for 7 ?t it remains constant (with slight 
favour to 7 = 1 /2). A large error for r/R > 1 for 7 = comes 
from the fact that the averaging area covers the entire object, the 
resulting normal is 0, and the information on the object geom- 
etry vanishes. If 7 ^ 0, the problem is avoided because normal 
vectors defined on facets on the boundary lying on the opposite 
side of the object have much lower weights. 



4. Wall shear stress in a channel flow 



In this section we test the accuracy of the WSS equation (|2]l 
combined with geometric normals in an inclined and bent chan- 
nel flows. 
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Figure 6: (color online) The maximum eiTor of geometrical and dynamical 
normal in a channel at various inclination angles. For geometrical scheme three 
r are tested. 




Figure 7: (color online) The maximum (a) and average (b) over 80 neighbour- 
ing sites of the error of WSS computed using various types of normal vectors 
in a channel at various inclination angles. Error bars in the average represent 
the standard en'or. 



4.1. Inclined channel 

A two-dimensional channel of an initial resolution 120 x 20 
lattice units (l.u.) was inclined at angle a to the x axis so that the 
radii of the inlet and outlet were kept constant. The analytical 
velocity profile (Poiseuille flow) was imposed on the inlet and 
outlet boundaries. The lattice viscosity was v = 0.01 and the 
maximum velocity in the channel center was Umax - 0.01 (lat- 
tice units). Using deviations from the analytical velocity field 
we found the steady state in less than 10"* iterations. We ran the 
flow for tan a e [0; 2] and calculated dynamical and geomet- 
rical normal vectors in a strip of length 80 l.u. (centered hori- 
zontally). In Fig.|6]the maximum angle between the numerical 
and exact vectors is given as a function of the channel incli- 
nation angle. The error of the geometric procedure decreases 
with r and is smaller than for the dynamic scheme. The most 
accurate results were obtained for tan a -Q and tan a - \. 

Next, we compute WSS using Eq. ^ and show that in an 
inclined channel WSS is essentially independent of the choice 
of the normal vector type (Fig. |7}. This implies that the largest 
contribution to the error comes from the velocity field near the 
staircase boundary. This fact is also visible in Fig.|6]where the 
large error of dynamic normals is caused by the velocity field 
inaccuracy. The maximum relative error of numerical WSS cal- 
culated at a single node is much larger than the error of the me an 
taken over 80 nodes. Therefore one could use an equation sim- 
ilar to Eq.m where a spatial average of WSS in a vicinity of the 
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Figure 8: (color online) Exact normals in a bent channel flow (short arrows) 
compared to dynamic (a) and geometric (b) normals. Boundary cells are empty 
and only pait of the channel is shown for both methods (the other parts are 
symmetric). 



wall element is computed: 

where w, are weights of the same form as used for normals. 
However, we leave applicability of Eq. (|5]l to more complex 
flows as an open problem. 

4.2. Bent channel 

Next, we investigate the flow in a two-dimensional bent chan- 
nel geometry defined by its inner and outer radii {R\ and R2, 
respectively). For the mesh of resolution W x // a site placed at 
(i, j) is of fluid type if R\ < |(/ - W/2, j)\ < R2, otherwise it is 
a no-slip site (or it is excluded). We take R\ - 10 and R2 - 20 
and the mesh resolution 45 x 25 l.u. and run the steady state 
simulation with v = 0.01 and analytical inlet/outlet boundary 
conditions. 

We calculated normals using both procedures and compare 
them visually against exact normals (Fig.[8]l. To understand the 
influence of the flow field near staircase wall on dynamic nor- 
mals the flow field is shown using streamlines (Fig. |8]a). We 
determined the angle between each individual pair of numeri- 
cal and exact vectors using the data from Fig. |8]and found the 
maximum error of the dynamic scheme to be around 25°. For 
the geometric scheme the error is kept below 10°. The average 
errors are calculated as 10.8° and 3.9° for dynamic and geomet- 
ric vectors, respectively. 

Then we computed WSS at the inner wall at various an- 
gular positions and plotted its components against exact solu- 
tions (Fig. |9]l. The numerical data follow the analytical solu- 
tions. Regions of larger deviations may be recognized (e.g. 
Tj, at a e (20,40)) mostly because WSS was calulated at the 
first node next to the boundary. An interesting observation may 
be done for the central region {a e (80, 100)) where introduc- 
tion of the geometric normals slightly decreased the error for 
Tx (similarly for Ty). The reason could be that the central part 
of the bent channel is horizontally flat in the grid representation 
(see Fig. [8]) and the velocity field follows the geometry. There- 
fore, the resulting dynamic normals are almost perpendicular 
to the surface whereas they are not in the original geometry (a 
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Figure 9: WSS components in the bent channel flow computed with Eq. J2) 
at various angular positions (2 l.u. from the inner channel wall). Results are 
obtained using geometric (squares) and dynamic (circles) normals. Empty and 
filled symbols represent the x and y component, respectively. Data are normal- 
ized by r,ef = max(T,). 



bent inner circle). This results in a large systematic deviation 
of the X component of WSS which is slightly decreased if more 
accurate geometric normals are used. 

5. A flow through the human abdominal aorta 

Finally, we test the procedures described in previous sections 
in a real situation and compute WSS in an abdominal aorta 
model based on patient-specific data. For this purpose we uti- 
lize data from the Virtual Family (VF) |51ll - a freely available 
library of high resolution anatomical human body models based 
on MRl measurements. We use Eq. Q to compute the geomet- 
ric normals and Eq. (|3) for WSS. 

5.7. Data preparation 

From several models provided by the VF we chose Duke - a 
34 year old man. Using the software included in VF we man- 
ually selected a part of his abdominal aorta 17.8 cm over and 
5.2 cm under the bifurcation point. We exported voxel data as 
.raw files that encode a tissue type with integers. Then we used 
the CVMLCPP library 115 211 ray marching algorithm to approx- 
imate the artery surface. The resulting mesh was of very poor 
quality, therefore further smoothing was done using MeshLab 
1 5311 . At this point we either voxelized the mesh to get various 
grid resolutions for the LBM or continued meshing until we ob- 
tained a three-dimensional unstructured volumetric mesh using 
NETGEN [541 . Next we applied the boundary conditions and 
exported the mesh in the neutral NETGEN format for further 
use in FVM software. 



5.2. Simulation 

First, we ran the steady state LBM simulation using the Sail- 
fish library adjusted for our purposes. We imposed a flat ve- 
locity profile on the inlet and zero pressure at both outlets 
of the model. Due to the tortuous geometry of the aorta the 
flat inlet velocity causes a macroscopic flow with a slightly 
lower effective flowrate. We used the following parameters: 
Mo,p = 0.0764 m/s (the velocity at inlet), /o,p = 0.03 m (the 
diameter of the inlet), v = 3 ■ 10"^ m^/s (kinematic viscosity 
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Figure 10: (color online) Wall shear stress calculated using LBM and FVM. The lighter the color, the larger WSS acts on the surface. 



of the fluid), Qp = 1000 kg/m^ (the fluid density). Thus, the 
characteristic time fo,p = lo,p/uo^p - 0.393 s and the Reynolds 
number Re=300 (based on an effective flowrate). The model 
was subdivided into 461 x 108 x 163 voxels (5mm resolu- 
tion). The number of grid nodes along the inlet diameter was 
A?=60, thus 6x = 0.0166. We chosed uih = 0.125, for which 
6t = uib ■ 6x = 0.00208 was found. The final lattice viscos- 
ity vib - 6t/6x~ ■ 1/Re - 0.0098 was found. We continued 
the simulation until the steady stade was reached (41 000 time 
steps) based on the observation of the velocity and pressure 
field convergence. The resulting distribution function that de- 
scribes the steady state was used in Eq. (|2]i to get WSS at 2 
l.u. away from the boundary. We used our geometric normals 
in Eq. (O as well as to move away from the boundary nodes. 
We first computed a dimensionless counterpart of T/b, namely 
Td - Qih ■ 6x~ I6f- ■ Tib - 63.7T/i. Then, the stress in [Pa] was 
obtained using Tp = Qp ■ lo.p/to,p ■ Td = SMTd- 

Second, we ran the same steady flow using finite volume 
Semi Implicit Method for Pressure Linked Equations (SIM- 
PLE) i55[| . We used simpleFoam solver from OpenFOAM Ii56i1 
with previously generated mesh built of 496 656 tetrahedral ele- 
ments. We applied a constant flowrate condition at the inlet with 
the flowrate Q=2.1206 ■ 10"'' [m^s] corresponding to Re - 300. 
The zero pressure and zero velocity gradient was applied at both 
outlets. We stopped the simulation after 1 second and confirmed 



that the convergence criteria was met (the initial and final resid- 
uals of pressure and velocity was smaller than 10"^). Then, we 
computed WSS using standard OpenFOAM tools. 

The resulting WSS in the abdominal aorta model computed 
with both methods is visualized as a color map (Fig.[T]and Fig. 
[To]). The resulting ranges of WSS were (6.74 ■ 10-^0.32) for 
LBM and (0, 0.42) for FVM. We found that the discrepancy 
in maximum values is due to the extreme values at the very 
boundary conditions, thus we encoded the WSS magnitude in 
the range of WSS E (0, 0.15) [Pa] only for better visualization. 
The profile of WSS (Fig.[TO]a) is taken along the white dashed 
line (see Fig. [10] b) from a gray scale projection of the results 
onto a plane. 

The peak value in the bifurcation area is around 0.12 [Pa] as 
observed in two branches behind the bifurcation point where 
the wavy profile of WSS may be easily recognized. The wavy 
profile of WSS seems to be mainly an effect of the geometry: 
the smaller profile diameter, the larger WSS. The quantitative 
comparison of WSS between LBM and WSS (Fig.[TO]a) shows 
an excellent agreement in the area not influenced by boundary 
conditions. 

Next, we compare the WSS magnitudes obtained for both 
methods with the experimental data measured using the laser 
photochromic dye tracer technique 115 711 . There, for similar con- 
ditions (the abdominal aorta of the healthy 35-year-old subject 



at rest and the steady flow) and the Reynolds number Re - 227 
(based on the inlet flowrate) the measurement values of WSS 
were in the range r e (0.06, 0.3) [Pa]. We found our numerical 
computations agree well with the experimental data. The re- 
maining diff'erence in the minimum value of WSS is most prob- 
ably an effect of the patient-specific geometry as we found the 
minimum of WSS in Fig. [TOh at y = 16cm, where the enlarge- 
ment of the cross-section of our aorta model is clearly notice- 
able. 

6. Discussion 

Comparison of normal vectors computed for three- 
dimensional objects with two diflTerent methods given in Sec. 
13.11 indicates an improvement in accuracy of our procedure 
against the dynamic scheme based on the velocity field. The 
improvement is significant e.g. close to the wall for a small lat- 
tice the error for dynamic normals is close to 90° but it remains 
below 10° in our method (Fig|6]l. Additional tests in the flow 
through an inclined (Fig.|6]l and bent (Fig.|8]i channel flows con- 
firm higher accuracy of geometric normals. Application of the 
geometric normals gives practically the same WSS as when its 
dynamic or exact counterparts are used in an inclined channel. 
In the bent channel, however, some differences are visible if 
separate stress components are analyzed (see Fig.|9). We found 
a region of an increased WSS accuracy based on the geometric 
normals (central part of the inner wall). We believe this obser- 
vations favorize the use of our scheme as it is simpler, more effi- 
cient (fewer algebraic operations) and provides the sense of the 
normal. We are aware that our method is nonlocal which may 
lead to complications in parallel implementation, this, however, 
might be solved e.g. by using algorithms that utilize fast shared 
memory of GPU processors. 

Our tests in a real patient-based geometry show an excellent 
agreement with a standard finite volume solver The WSS in the 
fluid flow through the abdominal aorta is both qualitatively and 
quantitatively the same in LBM and FVM (see Fig. [TOb . The 
slight difference between both methods at inlets and outlets is 
an effect of differences in the boundary conditions used. The 
data in Fig.[TO]a) give an impression of what is a typical length 
scale at which the imposed inlet boundary conditions play a role 
in hemodynamical simulations (e.g. we got around 10 cm from 
the inlet for the given setup). 

To get acceptable matching of the results in Fig.[TO]we moved 
the computation 2 nodes away from the wall (as suggested in 
M45i'l ). This procedure has, however, some limitations. First, 
while producing color-maps of WSS we noticed some missing 
nodes (see discontinuity in Fig. [TOb just above the bifurcation 
point). The reason for it is the following: if we start from two 
neighbouring sites and move some distance along their nor- 
mals, the destination nodes no longer have to be neighbours 
(especially if normals were pointing in different directions) and 
empty nodes at the surface may be visible. This artifact does 
not influence WSS values obtained. To produce decent visu- 
alization one could additionaly compute WSS for those empty 
sites, e.g. in the Fig. [1] we visualize WSS at a distances d=2 
and at d=0 which practically eliminated the problem of empty 



sites. Second, moving away from boundary may be unwanted 
in general because one could argue the quantity is no longer 
the wall shear stress but rather a near-the-wall shear stress. To 
make the use of the information at the wall surface only, we 
suggested averaging of WSS at the wall in the neighbourhood 
of a site for which WSS is calculated (see Eq. |5]i. Averaging 
decreased the error of WSS at almost each site of the inclined 
channel flow by more than 50% (see Fig.|7]i. To make the WSS 
averaging procedure a part of a complete protocol one should 
further test its behavior for smaller averaging radii and various 
three-dimensional geometries. 

We found that for a single model the preparation of FVM 
mesh described in Sec. 15. II takes around an hour without data 
segmentation (the Virtual Family implements this). In contrast, 
LBM could utilize the segmented MRI/CT data directly and 
skip the meshing process. We can imagine that, if the gray level 
was decoded correctly into some of physical variables used in 
LBM simulation (e.g. endothelial wall roughness) then these 
data could be used directly without segmentation. In this way 
one could turn the largest weakness of the LBM method (stair- 
case approximation of the boundary) into its strongest point in 
hemodynamical applications (no need for tedious data prepro- 
cessing), where our simple procedure for normals would be of 
high importance. 

7. Conclusions 

There are several reasons why one would want to calculate 
the wall orientation in LBM using the geometric procedure de- 
scribed here. First, our method gives a gain in accuracy over 
dynamic normals. Second, the computation overhead is much 
smaller as the averaging is done using a few simple algebraic 
operations rather than by solving an eigenvalue problem. Third, 
we showed that our procedure gives practically the same WSS 
as if exact normals were used. Fourth, the computation of geo- 
metric normals is - by definition - independent of the Reynolds 
number whereas dynamic normals are dependent on Re, which 
is unphysical. Fifth, it is suitable for the time dependent flows. 
Sixth, the given scheme provides the sense of normals that 
might be crucial in the calculations that require determination 
of separated wall shear stress components e.g. oscillating shear 
index and wall shear stress gradient. 
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Figure A. 1 1 : (color online) The overall stress T'"' and its tangential part t. 

Appendix A. The stress on arbitrary plane in LBM 

The stress tensor in fluids is a sum of pressure and viscous 
terms: 

o-ij^-pSij + a-'fy (A.l) 

The Cauchy formula gives the overall stress on the wall with a 
normal vector n (Fig. lA.llb : 



r. = O-ijilj. 



(A.2) 



The shear stress may be computed as the difference between the 
overall stress and its projection onto the normal: 



T; = (TijUj - {crkjnjnk)ni. 



(A.3) 



We insert Eq. dA.lb into Eq. (IA.3b . and since djktijrik - 1, the 
equation for WSS reads: 



T; = o-'fjiij - {o-[jnjnk)n 



(A.4) 



For the incompressible, Newtonian fluid the viscous stress cr'.. 



is a linear function of the rate of strain sif 



(A.5) 



Combination of Eq. ( IA.4I ) and Eq. ( IA.5I ) gives: 

T/ = 2^ [sijiij - {skjUjUkj n,) . (A.6) 

A definition of the rate of strain tensor in LBM, which uses the 
non-equilibrium part of the distribution function reads: 






2cjg 



^ai*-aj^ 



(A.7) 



where fl^"' - fa - fa' ^nd f"' is an equilibrium Maxwell- 
Boltzmann distribution. Finally, we combine Eqs. ( IA.6I I and 
( IA.7b and and get Eq. Q. 



Appendix B. WSS algorithm 



Ido 



-> dim - 1 do 

sumk + CakHink 



Algorithm 1 Calculation of WSS from Eq. (|2]i. 
dim <— 3 

C «- (yUw)/(cj£)) 

for / = — > ditn 
suma <— 
for all a do 

sumk <— 
for /t = - 

sumk 
end for 
sumj <— 
for j = ^ dim - 1 do 

sumj «— sumj + Cajtij (c„/ — sumk) 
end for 

Ja Ja Ja 

suma <— suma + fa"' ■ sumj 
end for 
T; <— C ■ suma 
end for 

^^55^ Irl 
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